library(haven)
II_Joo_Mukherjee_12 <- read_dta("II_Joo_Mukherjee_12.dta")
View(II_Joo_Mukherjee_12)
jm = II_Joo_Mukherjee_12
coef <- c(0.4455619,
0.0157448,
-0.3074662,
-0.0343053,
0.264737,
0.2530055,
0.2918376,
-0.1933579,
0.0040991,
-0.2345839,
-0.0237582,
-0.0040464,
-27.90933,
2.099829,
-0.0496096,
-0.9857114)
vcov <- as.matrix(read.csv("vcov.csv"))
lag <- 0
rebst <- 1
nego <- 0
cease <- 0
mob <- 1
lbd <- median(jm$logbd, na.rm=T)
lbd2 <- median(jm$logbd2, na.rm=T)
ind <- 0
terri <- 0
t <- median(jm$t1, na.rm=T)
t2 <- median(jm$t2, na.rm=T)
t3 <- median(jm$t3, na.rm=T)
censt.sim <- c(1,2,3)
mtime <- median(jm$dur_exist, na.rm=T)
sdtime <- sd(jm$dur_exist, na.rm=T)
set.seed(88110905)
B <- 3000
out <- matrix(NA, ncol=3, nrow=B)
ext <- matrix(NA, ncol=6, nrow=B)
for(i in 1:length(censt.sim)){
censt <- censt.sim[i]
pr.bs <- matrix(NA,B,1)
for(b in 1:B){
newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
bs.sum <- newcoef[1]*lag + newcoef[2]*censt.sim[i]*(mtime) + newcoef[3]*censt.sim[i] + newcoef[4]*(mtime) + newcoef[5]*nego +
newcoef[6]*cease + newcoef[7]*rebst + newcoef[8]*terri + newcoef[9]*mob + newcoef[10]*ind + newcoef[11]*lbd + newcoef[12]*lbd2 +
newcoef[13]*t + newcoef[14]*t2 + newcoef[15]*t3 + newcoef[16]
bs.sum2 <- newcoef[1]*lag + newcoef[2]*censt.sim[i]*(mtime+sdtime) + newcoef[3]*censt.sim[i] + newcoef[4]*(mtime+sdtime) + newcoef[5]*nego +
newcoef[6]*cease + newcoef[7]*rebst + newcoef[8]*terri + newcoef[9]*mob + newcoef[10]*ind +  newcoef[11]*lbd + newcoef[12]*lbd2 +
newcoef[13]*t + newcoef[14]*t2 + newcoef[15]*t3 + newcoef[16]
out[b,i] <- pnorm(bs.sum2) - pnorm(bs.sum)
if(i==1){
ext[b,1] <-pnorm(bs.sum2)
ext[b,2] <- pnorm(bs.sum)
} else{
if(i==2){
ext[b,3] <-pnorm(bs.sum2)
ext[b,4] <- pnorm(bs.sum)
} else{
if(i==3){
ext[b,5] <-pnorm(bs.sum2)
ext[b,6] <- pnorm(bs.sum)
}
}
}
}
rm(bs.sum)
rm(bs.sum2)
}
library(foreign)
library(stargazer)
library(ggplot2)
library(mvtnorm)
library(ggthemes)
library(reshape2)
library(margins)
library(survival)
library(survminer)
### load libraries ###
install.packages("margins")
library(margins)
for(i in 1:length(censt.sim)){
censt <- censt.sim[i]
pr.bs <- matrix(NA,B,1)
for(b in 1:B){
newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
bs.sum <- newcoef[1]*lag + newcoef[2]*censt.sim[i]*(mtime) + newcoef[3]*censt.sim[i] + newcoef[4]*(mtime) + newcoef[5]*nego +
newcoef[6]*cease + newcoef[7]*rebst + newcoef[8]*terri + newcoef[9]*mob + newcoef[10]*ind + newcoef[11]*lbd + newcoef[12]*lbd2 +
newcoef[13]*t + newcoef[14]*t2 + newcoef[15]*t3 + newcoef[16]
bs.sum2 <- newcoef[1]*lag + newcoef[2]*censt.sim[i]*(mtime+sdtime) + newcoef[3]*censt.sim[i] + newcoef[4]*(mtime+sdtime) + newcoef[5]*nego +
newcoef[6]*cease + newcoef[7]*rebst + newcoef[8]*terri + newcoef[9]*mob + newcoef[10]*ind +  newcoef[11]*lbd + newcoef[12]*lbd2 +
newcoef[13]*t + newcoef[14]*t2 + newcoef[15]*t3 + newcoef[16]
out[b,i] <- pnorm(bs.sum2) - pnorm(bs.sum)
if(i==1){
ext[b,1] <-pnorm(bs.sum2)
ext[b,2] <- pnorm(bs.sum)
} else{
if(i==2){
ext[b,3] <-pnorm(bs.sum2)
ext[b,4] <- pnorm(bs.sum)
} else{
if(i==3){
ext[b,5] <-pnorm(bs.sum2)
ext[b,6] <- pnorm(bs.sum)
}
}
}
}
rm(bs.sum)
rm(bs.sum2)
}
low <- as.matrix(sort(out[,1]))
moder <- as.matrix(sort(out[,2]))
high <- as.matrix(sort(out[,3]))
sort.out <- cbind(low, moder, high)
conf95 <- sort.out[(B*0.05+1):(B*0.95),]
pdf(paste(wd, "/Figure_files/Figure1.pdf", sep=""),width=8,height=8,paper='special')
boxplot(conf95, boxwex=0.5, names=c("Low", "Moderate", "High"), xlab="Central Command",
ylab="Change in Pr(Split)", cex.lab=1.5, cex.axis=1.5, cex=1.5)
abline(h=0, lty=2, col='red')
dev.off()
coef <- c(0.3732367,
0.0169961,
-0.355772,
-0.0414456,
0.2772929,
0.4089621,
0.3525784,
-0.1484045,
-0.0332967,
-37.63297,
3.684742,
-0.0841062,
-1.430015)
vcov <- as.matrix(read.csv("vcov_fig2.csv"))
lag <- 0
rebst <- 1
nego <- 0
cease <- 0
mob <- 1
ind <- 0
terri <- 0
t <- median(jm$t1, na.rm=T)
t2 <- median(jm$t2, na.rm=T)
t3 <- median(jm$t3, na.rm=T)
censt.sim <- c(1,2,3)
mtime <- median(jm$dur_exist, na.rm=T)
sdtime <- sd(jm$dur_exist, na.rm=T)
time.sim <- seq(min(jm$dur_exist), max(jm$dur_exist), by=0.5)
#############################
#### Figure 2a: High  #######
#############################
B <- 1000
set.seed(905)
out1 <- matrix(NA, ncol=3, nrow=length(time.sim))
for(i in 1:length(time.sim)){
time <- time.sim[i]
pr.bs <- matrix(NA,B,1)
for(b in 1:B){
newcoef <-  rmvnorm(1, mean=coef, sigma=vcov)
bs.sum <- newcoef[1]*lag + newcoef[2]*3*time + newcoef[3]*3 + newcoef[4]*time + newcoef[5]*nego +
newcoef[6]*cease + newcoef[7]*rebst + newcoef[8]*terri + newcoef[9]*mob + newcoef[10]*t + newcoef[11]*t2 + newcoef[12]*t3 + newcoef[13]
pr.bs[b,] <- pnorm(bs.sum)
}
out1[i,] <- c(mean(pr.bs, na.rm=TRUE), quantile(pr.bs, c(0.025, 0.975), na.rm=TRUE))
rm(bs.sum)
}
pdf(paste(wd, "/Figure_files/Figure2a.pdf", sep=""),width=10,height=8,paper='special')
par(mar=c(5,5,1,1))
plot(time.sim, out1[,1], type="n", xlab="Rebel Age", ylab="Pr(Split)",
ylim=c(-0.0005,0.3),
cex.lab=1.8, cex.axis=1.5)
lines(lowess(time.sim, out1[,1]), lty=1, lwd=2)
lines(lowess(time.sim, out1[,2]), lty=2, lwd=2)
lines(lowess(time.sim, out1[,3]), lty=2, lwd=2)
dev.off()
getwd()
